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» ' 

■ We generalize a proposal for detecting single phonon transitions in a single nanoelectromechanical 
' system (NEMS) to include the intrinsic anharmonicity of each mechanical oscillator. In this scheme 

two NEMS oscillators are coupled via a term quadratic in the amplitude of oscillation for each 

■ oscillator. One NEMS oscillator is driven and strongly damped and becomes a transducer for 
phonon number in the other measured oscillator. We derive the conditions for this measurement 
scheme to be quantum limited and find a condition on the size of the anharmonicity. We also derive 
the relation between the phase diffusion back-action noise due to number measurement and the 
localization time for the measured system to enter a phonon number eigenstate. We relate both 
these time scales to the strength of the measured signal, which is an induced current proportional 
to the position of the readout oscillator. 
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^3 ' I. INTRODUCTION 



With device fabrication in the submicron or nanometer regime, it is possible to fabricate mechanical oscillators with 
very high fundamental frequencies and high mechanical quality factors. In the regime when the individual mechanical 
quanta are of the order of or greater than the thermal energy, quantum effects become important. Recently, a high- 
er frequency mechanical resonator beam that operates at GHz frequencies has been reported 1] . Unlike quantum optical 
systems where extremely high frequency oscillators, vacuum environments, zero temperature, and well-isolated systems 
are the usual setup, solid state systems normally exist at finite temperatures and interact with their surroundings. 
For a resonator operating at the fundamental frequency of GHz and at a temperature of lOOmK, on average only 20 
vibrational quanta are present in the fundamental mode. An interesting question is: can we observe quantum jumps, 
i.e., discrete (Fock or number state) transitions in such a true mechanical oscillator in a mesoscopic solid system 0, 
as the mechanical oscillator exchanges quanta with the outside world or environment? In order to observe quantum 
jumps, one needs to design a scheme to measure the phonon number of the oscillator so that the oscillator will stay 
in a certain phonon number state long enough before it jumps to another phonon number state due to the inevitable 
interaction with its environment, usually through linear coupling to the oscillator position. 

To achieve a quantum mechanical phonon number measurement of a mechanical oscillator, conventional measure- 
ment methods, such as the direct displacement measurement, [j| cannot be simply applied since the observable (i.e., 
the number of phonons in the oscillator) does not commute with, for example, the position or displacement operator. 
Thus, naively attaching a readout transducer to the mechanical oscillator results in inaccurate subsequent measure- 
ments due to back action. One thus must make sure that the transducer that couples to the mechanical resonator 
measures only the mean-square position, without coupling linearly to the resonator's position itself 0- 

Some preliminary experiments in this direction have been conducted |22|. They use a second, driven mechanical 
oscillator (oscillator 1 in Fig. ^) as the transducer to measure the mean-square position of the system oscillator 
(oscillator in Fig.0. Hereafter, we use the notations of the "system oscillator" and "ancilla oscillator" in the text, 
but keep and 1 as subscripts in the mathematical notations. The basic idea is that the non-linear, quadratic-in- 
position coupling between the two oscillators shifts the resonance frequency of the ancilla oscillator by an amount 
proportional to the phonon number or energy excitation of the system oscillator. This frequency shift may be detected 
as a phase shift of the oscillations of the ancilla oscillator with respect to the driving, when driven at a fixed frequency 
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near resonance. Also, the ancilla oscillator needs to have sufficient sensitivity to resolve an individual quantum jump. 

In the analysis of this measurement scheme presented by Santamore, Doherty, and Cross |£g, self-anharmonic terms 
xj, in the two mechanical oscillators were neglected due to the smallness of the coupling coefficients compared to their 
harmonic oscillation frequencies, where Xi, i = 0, 1 is the displacement of the oscillators position from equilibrium. 
Since the self-anharmonic terms are of the same order as the nonlinear coupling term Xqx\, it is important to include 
those terms and analyze the effects on the proposed measurement scheme. 

In this paper, we extend the work of Ref. 5] and investigate the effects of self-anharmonic terms on a phonon number 
measurement. Due to the higher order self-anharmonic terms, the adiabatic elimination method used in Ref. 5] may 
not be straightforwardly applied even with the assumption of a heavily damped ancilla oscillator due to measurement. 
Here we take a slightly different approach. As the ancilla is assumed to be heavily damped, it will relax very rapidly 
to its steady state within a timescale on the order of typical response time of the system oscillator, and will appear 
to the system oscillator effectively as a "bath". To see the consequences of a rapidly decaying ancilla oscillator on 
the dynamics of the system oscillator, we use the quantum open systems approach to find the master equation for 
the reduced density matrix of the system oscillator. In obtaining the master equation, the correlation functions of 
the "effective bath" (or the ancilla oscillator) are calculated using the generalized P-representation approach [f| • The 
generalized P-representation approach has the advantage of removing some of the unnecessary restrictions imposed 
in Ref. 0. 

We find that in the presence of self-anharmonic term, xf, of the ancilla oscillator, the effect of increasing driving 
strength and self-nonlinearity tends to shift the resonance frequency, increase the peak value and decrease the width 
of the response of the peak of (T/To) (see Fig. [2J. The quantity (T/Tq) is the ratio of the back-action diffusion 
coefficient (or decoherence rate) T [see Eq. Q55[l] and its value T at zero self-anharmonicity and zero detuning. If 
the damping of the ancilla oscillator is much larger than the effect of the self-anharmonic term, the overall effect of 
self-anharmonic term on the phonon number measurement is small. Finally, we show that the induced electromotive 
readout current 8] from the ancilla oscillator provides information on the phonon number of the system, even in the 
presence of higher order anharmonic terms, and we obtain the relation between the current and the measured system 
observable. 

In the next section, we discuss briefly the measurement scheme and Hamiltonian, and obtain the master equation 
for the model described above while keeping higher order self-anharmonic terms. It turns out that the master equation 
we obtain requires two-time correlations of the ancilla oscillator operators. Section ll VI deals with this issue. We find 
one-time and two-time correlation functions of the ancilla. In Sec. [V] we examine the effect of the self-anharmonic 
terms on the dynamics of the system oscillator from the master equation of its reduced density matrix. In Sec. I VII we 
obtain the dependence of the measurement current on the measured system oscillator observable, the phonon number. 



II. HAMILTONIAN AND THE MASTER EQUATION 



A. Proposed scheme 



thermal noise 



thermal noise 



1 



mechanical 
oscillator 



T 



A. 01 Xg Xj 



1 



mechanical 
oscillator 1 



T 



driving 



readout 



damping 



damping 



FIG. 1: Schematic of phonon number measurement for a coupled mechanical oscillator. Oscillators and 1 are anharmonically 
coupled with coupling strength Aoi. Both oscillators are subjected to thermal noise injection and dissipation. The oscillator 1 
is driven and a readout apparatus is attached to it. 
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Our model consists of two mesoscopic scale mechanical bridges with rectangular cross section. One serves as a 
system oscillator (oscillator in Fig. to be measured. The other is used as an ancilla oscillator (oscillator 1 in 
Fig. nj, and is part of the measuring apparatus. Details of the scheme have been already discussed in Ref. p|. 
A schematic illustration is reproduced in Fig. ^ These mesoscopic-size elastic bridges or beams with rectangular 
cross-section are connected by a device that transmits only one of the flexing modes of the system oscillator to the 
ancilla oscillator. As a result, these two resonators are anharmonically and symmetrically coupled (for experimental 
progress toward the scheme, see Ref. |1J). We label the measured system oscillator with subscript and the ancilla 
oscillator with subscript 1, with corresponding resonant frequencies of the two flexing modes labelled as luq and u>i, 
respectively. The ancilla oscillator is driven at frequency u>d with strength e(t). A measuring apparatus is attached 
to the ancilla oscillator. The whole structure is subjected to the thermal bath environment. The interaction of the 
system oscillator with the thermal bath causes thermal dissipation and excitation of the system oscillator, which 
results in random-in-time transitions between phonon number eigenstates (i.e., quantum jumps). A change in the 
energy of the system oscillator appears to the ancilla oscillator as a shift of the resonant frequency via the anharmonic 
coupling. This frequency shift may be detected as a phase shift of the oscillations of the ancilla oscillator with respect 
to the driving, when driven at a fixed frequency near resonance. 



B. Model Hamiltonian 



The free Hamiltonian for the two bridge oscillators and 1 is 

Hft ee = fkijQa'a + fuvib' i b, (1) 

where a^,a are creation and annihilation operators for oscillator 0, respectively, and similarly, b\b for oscillator 1. 
The ancilla oscillator is driven at frequency Wd with driving strength e, 

#drive = he cos [uj d t] (6 f + b) , (2) 

In the interaction picture, the driving term becomes 

Hl rive = 2he(tfe +iS » + be- iS »), (3) 

where 5u> is detuning between the ancilla resonant frequency and the driving frequency, lo\ — U)d- 

The two oscillators are coupled anharmonically through the special coupling device that controls and allows only 
one type of strain (the longitudinal stretch) to pass to other oscillator. Beyond the linear elasticity theory, the two 
flexing modes, which are perpendicular to each other, are coupled. Expansion of the elastic energy with respect to 
the strain tensor is taken up to second order. The next term, cubic in the elastic energy, gives quadratic terms in the 
equation of motion [9L Il0|. Since the coupling of the two modes of the two beams is symmetric, and since the two 
modes are not coupled at the linear level, the first order in coupling is xf s x\, where Xi is the displacement operator. 
So we expand the anharmonic terms up to first order in coupling and obtain 

Hanh = h (a £o + -Vra^o + hxl + Xnxfj , (4) 

Vanh = HXoixlxf, (5) 

where A^ is the coupling coefficient. The high frequencies of the resonators, i.e., (lo — u>i) much larger than A01 and 
their damping rates, allows us to use the rotating wave approximation. Thus we write the anharmonic terms as 

Hanh = ftA o (a 1 a) 2 + HX U (tfb) 2 , (6) 
Vanh = hXaiaJatfb, (7) 

where we have defined the standard raising and lowering operators for the oscillators: a — \J mou>o/2h xo + 
i^/l/2hmoLiJQ po, o-^ is the Hermitian conjugate of a and similarly for b and b^ with the subscript replaced by 
1. We have also introduced new coefficients A's (without tildes) which all have the same dimension of frequency. 

The coupling term Xoia^ab^b commutes with the observable a^a, enabling a quantum non-demolition (QND) mea- 
surement. The terms Xoa'a and Xib^b shift the resonance frequency by a constant amount, so we have absorbed these 
quantities into ujq and uj\. The terms fi.Aoo(a^a) 2 and hXn(b'b) 2 are analogous to Kerr non-linearities in nonlinear 
optics. Since these terms commute with the measured observable (a* a), they will not change the system phonon 
number eigenstates; however, the Kerr effect causes an intensity dependent phase shift. Unlike a coherent state, in 
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which this effect results in rotational shearing, a thermal state will not be affected by phase shift, due to its rotational 
invariance. 

As for detecting phonon number in the system oscillator, we adapt a magnetomotive detection scheme suggested 
by Yurke et al |JQ ^M- The voltage developed is proportional to dx\/dt, where x\ is the displacement of the 
beam from its equilibrium position. The current induced by this voltage is monitored by phase lock-in amplifier. An 
experimenter monitors the amplitude of the current and its phase with respect to the driving current that is set to a 
frequency near resonance. The details of the relation between the measured current and the phonon number of the 
system oscillator are derived in Ref. (|5|). 

There are two physically distinct environments in the model: the thermo-mechanical environment of each oscillator 
and the electronic noise environment of the electrical system that ultimately provides information on the motion 
of the ancilla. The environments are modelled as thermal baths, each consisting of an infinite number of harmonic 
oscillators. The couplings between the oscillators and the thermal baths are considered as weak, linear and Markovian; 
thus we use the rotating wave approximation. The Hamiltonian of the baths and their coupling to the oscillators can 
then be written as: 

Hbath = u s,iBl tn B s>n , (8) 

s n 

Bm ) , (9) 

where s runs over three different baths: the thermal baths coupled to the system oscillator (BO) and ancilla oscillator 
(Bl), and the electronic (measurement) bath coupled to ancilla oscillator (Bm). The operator 

O s = ^ g s B s<n (10) 

n 

consists of bath operators, and the coupling to the bath modes is given by the coefficients g s (u> n ). 



C. Master equation 



Using the standard technique for open quantum systems, we first obtain the master equation for the joint density 
matrix of the two oscillators, R, by tracing out the bath variables: 

— — = — iuJa [a^a,i?] — iXoo (a 1 a) ,R 



5u [6+6, R] - ie [tf + b, R] - iXu (tfbf ,R - i\ 01 [aW&, R] 

> (N + 1) V [a] R + vN V [a f ] R + K (Ni + 1) V [b] R + nN{D [6 f ] R, 



where 



V[0]R = 20i?O t - (0 ] OR + RO^O), 
V [O f ] R = 20^RO - (OO t R + ROO r ). 

are defined for arbitrary operators O and R. The damping rate of the system oscillator v is given by 

v = itqbo (uj q ) \g B o (ojo)\ 2 ■ 



(11) 

(12) 
(13) 

(14) 



It is related to the quality factor Qq of the system oscillator by v = cjo/2Qo- We have combined the damping rates /i 
and r], due respectively to thermal bath and measurement on ancilla oscillator, into k = /j + 77, where 



/1 = irg Bm (ujx) \g Bm (ll>i)\ 2 , 
?7 = ttq b1 (u>i) \g B1 (wi)| 2 . 

Here g s (w) is the density of states of bath s at frequency u>. The Ni are the Bose-Einstein factors: 

1 



(15) 
(16) 



(17) 
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and Ni = (tjNi + fiN m ) /k, where 



Nt = -— !: , N m = ^ , (18) 



with /3; = (fceTs) 1 and T s the temperature of bath s. In Eq. 1|11|) . the first and the second lines are the free 
Hamiltonian and non-linear Kerr effect terms for system and ancilla oscillators, respectively. The third line in Eq. 
(|ll|l is associated with the anharmonic coupling, and the last two lines are consequences of the interactions with 
thermal baths. 



III. EFFECT OF HEAVILY DAMPED ANCILLA OSCILLATOR 



To proceed further towards a master equation for the reduced density matrix for the system oscillator alone, the 
ancilla oscillator is assumed to be heavily damped due to measurements, i.e., k 3> Xij,v. In this case, the ancilla 
oscillator will relax very rapidly to its steady state and appear to the system oscillator as a "bath". In fact, if 
An <C lui and Aoi <C n, the ancilla oscillator in Eq. (|llf> will remain near a thermal steady state with average number 
N\. However, we will relax the condition An <C u>i and treat the interaction Aoi term pertubatively. 

To see the consequences of the rapid decay of the ancilla oscillator on the dynamics of the system oscillator, we use 
perturbation theory and expand the interaction Hamiltonian Hj (t) = Xoia^ab^b up to second order, and trace out the 
ancilla oscillator variables. This implies that we need to calculate the relevant steady state averages and correlation 
functions for the ancilla oscillator in the presence of the anharmonic term An (b^b) . 

In this case, the master equation for the reduced density matrix p(t) for the system oscillator alone can be written 

as 

= -iu [at a, p{t)] - i\ oa [(a 1 ' af , p{t)\ +v(N + l)T> [a] p(t) + vN^V [a f ] p(t) 

- ill! [Hj (t) , R cS (t)} - f Tn [Hj (t) , [Hj (t') , R cS (t)}} dt' , (19) 

Jo 

where R e g(t) w p{t)p\(t) is the effective joint density matrix of the two oscillators under the approximation that the 
ancilla oscillator is heavily damped, and pi(oo) is the steady state density matrix operator for the ancilla oscillator. 
Explicitly, the second term of the last line of Eq. (|19|l can be written as 

t 

Tn [Hj(t),[Hj(t'),R cS (t)}}dt' 

- (Aoi) 2 / a) a (t) a) a (f) p(t) (b ] b (t) b ] b {t')) dt' + (A i) 2 / a)a (t) p{t)<Ja (*') (6 f 6 (f) tfb (t)) dt' 



+ (A i) 2 / a)a{t')p{t)a)a{t)(tfb{t)tfb{t'))dt' - (A i) 2 / p{t)a ] a{t') a) a{i) (b ] b{t')b ] b {t)) dt' . (20) 



The exact correlation functions of the ancilla oscillator are not easy to evaluate because of the presence of the 
anharmonicity, the driving, and the decay terms. However, one can make an expansion of the state of the ancilla 
oscillator around its steady state and linearize the fluctuations, assuming them to be small 0,^3- 

Define the steady-state mean field amplitudes as (b) = (3q. The operator b can be written in terms of small 
fluctuations about the steady state mean value as 

b (t) =Po + h (t) . (21) 

Then, keeping terms up to quadratic order in bi,b\, the interaction Hamiltonian Hj = Aoia^ ab^b becomes 



Hj = Aoicra 



\Po\ 2 +P* bi(t) + P b\(t)+b\(t)b 1 (t) . (22) 



The first term in Eq. (122H contributes to a shift in the resonant frequency of the system oscillator by a constant 
amount and can be combined with the free Hamiltonian. Inserting this expression back into the first term of the last 
line of Eq. (|19|) gives the first order expansion term 

-iAoilh [ot a (t) tfb (t) , Reg®] = -i\ 01 [a^a (t) , p] , (23) 
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where we have used the fact that averages of fluctuation fields vanish, i.e., 

(h) = {b\)=0. (24) 



Now we turn our attention to the second order term, Eq. (|2U[1 . Note that since k ^> v, the phonon number a^a^) 
of the system oscillator changes with time on a time scale much larger than b'b(t) of the ancilla oscillator. So we can 
approximate (va(t') ~ a/a{t) in Eq. I|20(l and pull the system oscillator terms outside of the integral. Then Eq. I|20() 
becomes 

t 

Tn .[^(O.iM*)]]*' 







~ (^oi) 2 \a)a (t)p(t)a*a(t)- (a^a{t)) 2 p(t)} J (B (t, t')) dt' 

+ (A i) 2 {at a {t)p{t)a) a {t) - p(t) (a? a (i)) 2 } f (B (t' , <)) dt' (25) 



where 

(B (t, 0) = Wo f (h (t) h (t')) + |/3 | 2 (fti (t) 6l (i')) + |/3 | 2 (&I (t) 6i (*')) + (A,) 2 (&I (t) &I (*')) (26) 
(B (f , *)) = (/3 *) 2 (6x (f ) 6i (t)) + |/3 | 2 (&i (*') b\ (t)) + |/3 | 2 (&+ (*') 6 X (t)) + (A,) 2 (&J (f ) &1 (*)) (27) 

and higher order fluctuation terms than b\ are ignored. The linearization transforms the second-order correlation 
functions of the ancilla operators, (b^ b(t)tf b(t')) and (b^b (t') b^b (i)), into first order correlation functions of fluctuation 
fields: (b\ (t) h (t')), (b\ (t) b\ (f )) (b\ (t) h (f )), and (&J (f) b \ (*')>. 

IV. ONE-TIME AND TWO-TIME CORRELATION FUNCTIONS OF ANCILLA 

In this section we calculate the one-time and two-time correlation functions of the ancilla oscillator. For this 
purpose, first we need to calculate the one-time correlation functions of a single driven anharmonic oscillator. We will 
follow the method of Drummond and Walls |14| . who obtained one-time correlation functions. Then we extend their 
method to calculate two-time correlation functions. 

The master equation for the driven, anharmonic ancilla oscillator interacting with the thermal bath is given by 

= -iSu [tfb, Pl {t)} - te [ftt + b , Pl (t)] - iX u [(6+6) 2 , Pl (t) 
+ k(N 1 + l)V[b] Pl (t) + KNiV [6 f ] pi{t). (28) 

where p\ is the density matrix of the ancilla oscillator and Su> = u)\ — u>d is the detuning, with u>d the driving frequency. 
The exact steady-state one-time correlation functions for a system with master equation Eq. (|28l) at zero temperature 
were given in Refs. |?l Il4|. in a discussion of optical bistability of a coherently driven dispersive cavity with a cubic 
nonlinearity in the polarizability of the internal medium. At finite temperature, no exact solution has been found. 

Our first objective is to derive a stochastic differential equation from the quantum master equation. Representing 
a density matrix in a coherent state basis is useful in systems described by Bose operators b\ 6. Due to the presence 
of the non-linear, self-anharmonic term, we will use the generalized P-representation introduced by Drummond and 
Gardiner @ to preserve the positivity of the Hermitian density operator. 

Using the above transformations, the Fokker- Planck equation corresponding to the master equation Eq. I)28|l can 
now be written as 

dt P (?) = { 80 [ (K + l5uJ + * All) & ~ 2lXlll32a + l€ \ ~ lXl1 d^P 2 

+ fa [(« - ife - iA ") a - ^n" 2 / 3 - ie ] ~ lXn ]^ a2 + 2KNi dpdc:\ p (p) ■ (29) 

The argument of the generalized P function is ft = (/3, a) T . The correspondence principle between operators and 
c- numbers is as follows: j3 <-> b and a <-> b^ . However, ((3, a) are not complex conjugates. Drummond and Gardiner 
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have shown [6( that the Fokker-Planck equation in /3 can be transformed to a stochastic differential equation with 
positive definite diffusion [2^ . They found that the stochastic differential equations in the Ito calculus corresponding 
to Eq. (HJ are 



(30) 



where £i and £2 are random Gaussian functions, so that (3 and a are complex conjugate in the mean j24j. This 
stochastic differential equation is non-linear and not solvable as it is. However, it is reasonable to use a small noise 
expansion and linearize the fluctuations about the steady state of the mean field amplitudes. Thus we write (3 in 
terms of the mean amplitude and first order expansion of the fluctuation, 



d 


' P~ 




—ie — f3 (k + iSui + iXn + 2iXii/3a) 




-2iX n f3 2 


2kN 1 


dt 


a 




it — a (k — iSuj — iXn — 2iAna/3) 


+ 


2kN x 


2i\na 2 



- 1/2 


'6" 







/?(*)=#, + A (t) 

where /?o is the steady-state mean amplitude of /3 and is given by 



i{5oj + X n + 2X 11 \f3 \ 2 ) + k 



(31) 



(32) 



and Pi is the zero mean fluctuation amplitude. We have a similar expression for a. Thus (3q and «o are complex 
conjugate to each other (i.e., Poao — \ao\ 2 = \/3q\ 2 = tiq)- Then to first order in the fluctuations, the fluctuation 
amplitude vector /3i = (/3i,ai) T obeys a stochastic differential equation 

d 



Of 



h (i) = -A-& (0 + D 1 / 2 (0o)i(t), 



(33) 



where £ = (£i,£2) T is the noise vector, A is the linearized drift matrix and D is the diffusion matrix evaluated at 
[3 = $q. The matrices A and D are 



k + i6u> + iXn + 4zAii?io 2zAh/?q 

— 2iAii«Q k — iSu — iXn — 4«Anno 



and 



D 



-2iXuf3 2 2kN x 
2kN x 2iXual 



(34) 



(35) 



The one-time correlation matrix can be calculated using the method of Chaturvedi, et al fft. ITl. HHITcl : 

r /«2\ /«,„„\i 

c (*,*) = 



(Pi) (Piai) 
(ai/3i) (a?) 



1 

A 2 " 



4iAnn ) (2iVi 

|2 



where 



-iXuP 2 (k — iSu — iAi 
iVi \k + iSuj + iXn + 4iAnno| 2 + 2Af 1 ng 

+ A 2 , 



1 



Ni \k + iSuj + iXn + 4zAnn | 2 + 2A 2 i n 2 , 
zAuOq (k + iSuj + iXn + 4iAnn ) (2Ni + 1) 



^2 _ ,„2 1 A 2 



A 2 = (5lu + An) 2 + 8 [Su + An) Aun + m^n 2 ,. 
We now derive an expression for the two-time steady state correlation matrix 



C(t,f) = 



(A(i)A(i')) C9i(t)ai(f)) 
(ai (t)/?i(*0) (01 (i) 01(0) 



For t > t' 



C (f, t') = exp (-A (t - t')) C (t, t) 



(36) 



(37) 
(38) 



(39) 



(40) 



and for t < t' 



C (t, t') = C (t, t) exp (-A T (t' - *)) 



(41) 
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Let us define M(i,i') = exp[— A(t — t')]. The matrix M can be calculated as follows. Let the matrix U = (ui,u 2 ) 
diagonalize A with eigenvalues A±. The eigenvalues for this 2x2 matrix can be found from the characteristic equation: 



Tr (A) ± J [Tr (A)f - 4 dot (A) 



k ± iAi, 



(42) 



We then obtain the matrix M as 

~exp(-A + (t-t')) 



M(t,t') = U 









exp(-A_ (t-t')) 



U 



1 

2AT 



(A x - c) e-^-C*-*') + (Ax + c) e- A +(*-*') 2A n /3 2 [-e-M'-*') + e - A +(*-'') 
2A u og \e~ x - (*-*') - e -^+(*-*')l (A x + c) e~ A - ('-*') + (A x - c) e - A +(*-*') 



(43) 



where c = 4An?in + <5cj + A xl . The two-time correlation matrix Eq. I|39[). then follows directly from Eqs. I|4(J|) (|41ll and 
(gSJ, as well as the fact that exp[-A(f - t')} = M(f, t') and exp[-A T (t' - t)} = M T (t', t): 

The detailed expressions of the two-time correlation functions are shown in the Appendix. We note that in the 
P-representation, the c-number time correlation function corresponds to a normally ordered time correlation function 
of the operators; thus the correlations above do not correspond to all the two-time correlation functions we need to 
find. For non-normally ordered time correlation functions, some care needs to be exercised. Using the procedure 
described, for example, in Refs. |l6l.ll7j. we obtain the following operator to c-number correspondence: 



(h (t)h (0) 


= (p 1 (0), 




(44) 


h (t)b\(t>)) 


= (/?! (t)ai (0) ■+ 


M n (t,t'), 


(45) 


b\(t)h (0) 


= <ai (t)j9i (t')), 




(46) 


b\(t)b\(t')) 


= (ai(t)ai(t')) 


-M 21 (t,t'), 


(47) 


(bi (t')h (t)) 


= (/3 x (t)# (t')> + 


M 12 (t,t'), 


(48) 


h (t')b\(t)) 


= (ai (t)/3i (0)4 


M 22 (t,t'), 


(49) 


b\(t')h (<)) 


= (/?! (t)oi (0), 




(50) 


b\(t')b\(t)\ 


= (ai (t) ai (*')> , 




(51) 



where Mij(t,t') are the matrix elements of the matrix M(t, t'), Eq. 14311 . 



V. MASTER EQUATION FOR A REDUCED DENSITY MATRIX 



Having found the one-time and two-time correlation functions, we can now evaluate Eqs. (|23|l and (|25|l and obtain 
the master equation for the reduced density matrix of the system oscillator as: 



where 



^ = -i (wo + A) [aVp] - iQ \\a'a) . r 
+ v (Nq + 1)T> [a] p + vNqD [a)] p 



— r \a'a, la'a, p\] 



Aqi 



n Q + — (^Ni \k + iSuj + i\n + 4iAnn | 2 + 2Ai 1 n 2 ] 



6 = Ann + 



A 2 



(Sui + i\n + 2Ann ) , 



-j^-Kn (2Ni + 1) \k + iSuj + iXu + 4iAiin | 2 — 4Ann (5w + A n + 3Ann ) 
^(2N 1 + 1). 



(52) 

(53) 
(54) 

(55) 
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We have set no = |/?o| 2 > aud A 2 is defined in Eqs. i|37|) and l|38(l . In obtaining the last line of Eq. I|55[l . we have used 
Eq. (22). 

In Eq. I)52[l. A in the first term is the resonant frequency shift due to interactions. The second term is the Kerr 
non-linear phase shift, with coefficient <d depending on the anharmonicity of both oscillators Aoo and An, as well as 
the detuning of the ancilla oscillator. The parameter T is the phase diffusion coefficient or decoherence rate, associated 
with back-action due to an effective measurement of a^a. Physically, due to monitoring, the system would localize or 
collapse into a phonon number eigenstate on a time scale of order r . The measurement time that is needed for the 
measurement apparatus to distinguish one state from the next is also proportional to r . The last two terms in Eq. 
(|52|l . can be derived from the thermal coupling to the system and are responsible for the quantum jumps. In the case 
when v = 0, the conditional master equation of Eq. (|52|) will describe a QND measurement of the system oscillator 
phonon number. The time the system stays in a given phonon number state before making a transition due to either 
excitation or relaxation is proportional to v~ x . To be a good quantum measurement of a phonon number state, we 
want the system's dwelling time to be long compared to the time necessary to determine which number state the 
system is in, i.e., (Y/v) ^> 1. 



A. Effects of the anharmonic terms 



From Eq. (|52|l . we notice several important points. Firstly, in the case of no detuning and no non-linear self- 
anharmonic terms (i.e., Slu = 0, Aoo = An =0), we have 

A = A 01 [Nx + (e/n) 2 ] , (56) 
6 = 0, (57) 

r= A °^ 2 f + (58) 

These results agree with the results of a simpler model discussed in Ref. Q, using a slightly different adiabatic 
elimination approach. 

Secondly, the steady state solution Eq. (|32J) of Eq. H^Qjl gives 



|e| 2 = n Q 



k 2 + (Slu + An + 2An?T.o)" ■ (59) 



Equation l|59|l has an analogy to a classical anharmonic oscillator Bistability due to a Kerr nonlinearity is a 
well known phenomenon. Classically the oscillator will take one or the other of the stable solutions. Using Hurwitz 
stability criterion, to obtain stable solution for Eqs. (f 33|> — l|35|) it is necessary to have 

Tr(A) > 0, (60) 
Dct(A) > 0. (61) 

For the matrix , Eq. J23J), gives Tr(A) = 2k > for a dissipative or loss mechanism. Therefore the threshold points 
are determined by Det(A) = A 2 = 0. However, in the quantum regime at zero temperature, bistability appears only 
during transient period and does not exist in the steady state. [6|,|3 We, nevertheless, note that the linear theory that 
we use to calculate the steady state correlation functions at finite temperatures would break down at the instability 
points. 

Secondly, from Eqs. I|52l) and (|55|l . we see that when Slu = 0, the condition k ^> An makes the effect of the non-linear 
self- anharmonic terms in A and T very small, which justifies the assumption of neglecting An in Rcf. 5]. However, 
our calculation allows us to do a quantitative analysis without making this assumption. 

The value of the phase diffusion coefficient T (as compared to the damping rate v) is important to the phonon 
number measurement scheme and to the observation of quantum jumps. To see the effects of self-anharmonicity, 
driving and detuning on the phase diffusion coefficient T compared to its value To at zero self-anharmonic coupling 
and zero detuning (An = e = 0) @, we plot their ratio 

£ - £ (62) 

in Fig. |21 Note that T diverges at A 2 = 0, which are the instability points where the linear theory is not valid. The 
parameters (in units of k) in Fig. [21 are chosen so that the ancilla oscillator is away from these points. For example, 
if we were to increase further the driving strength in the dot-dashed line plot of Fig. [21 to e = 1.2, say, the ancilla 
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FIG. 2: The ratio of (T/To) as a function of detuning at different values of driving strength and self-anharmonicity (Kerr effect 
coupling). The parameters are presented in unit of the damping rate k. 



oscillator would then be in the instability regime. When the nonlinearity An is small, the solid line plot in Fig. [21 
shows the linear resonance of small driving. The dotted, dashed, and dot-dashed line plots illustrate that increasing 
the driving strength and the nonlinearity tends to shift the resonance frequency, increase the peak value and decrease 
the width of the peak of (r/ro). 

Carr and Wybourne have estimated an anharmonic coefficient Xa for a beam with rectangular cross-section^^: 

A, = ^- 22 f 5 , (63) 

where B is the bulk modulus, pi is the mass density, L, w,t are the dimensions of the beam: length, width, thickness, 
respectively. A simple estimation of k and An using realistic values for a mesoscopic mechanical oscillator reveals 
that An is many order of magnitude smaller than n. 



VI. MEASUREMENT CURRENT 



In the measurement scheme, we do not observe the phonon number of the system oscillator directly. Rather we 
perform a phase sensitive, 'homodyne', measurement on the quadrature (b + tf) of the ancilla oscillator. It is therefore 
important to show that an observation of the average current (7) = ^/2]l(b+b ' ) indeed corresponds to a phonon number 
measurement of the system oscillator. We anticipate that the average measured current of the ancilla oscillator is 
proportional to the average phonon number in the measured system oscillator. In addition we need to show that the 
coefficient of proportionality is related to the localization rate, which determines how long it takes to distinguish one 
number state from the next. Thus a strong signal corresponds to a rapid localization rate. Furthermore we expect 
that the localization rate is proportional to the backaction induced phase diffusion coefficient T, so that the better 
the measurement, the larger is the back action noise. 

To demonstrate this, firstly we use the Hamiltonian to obtain the quantum Langevin equation for the ancilla 
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oscillator operator b: 



— = — ie — iSujb — iX^\a) ab — i2X\\b^bb — rjb (t) — y/2rjB^ n (t) 



dtf . r 

— = ie + iSutf +«A ia t a6 t + i2An& t fe t fe- rjtf (t) 



2^1 (*) 



lib (t) - ^2nD in (t) 
/i&t ( t ) _ (t) 



(64) 
(65) 



where £?i„ (t) is the input noise 17] . The steady state [(d&/<i£) = 0] average of (b) — /3q for the ancilla oscillator in 
isolation (i.e., with Aoi = 0) is given by the same expression as Eq. (|32[1 . Linearizing around the steady state, renaming 
the operator describing the quantum fluctuation as b\(t), and assuming that a* a do not to change appreciably over 
the typical time scale of the ancilla oscillator, we obtain 



dbx 
~dt 
db[ 
~~dt 



~i(5uj + An)bi + 2An(2ri 6i + /3$b\) + Xoifoala - nb x 
i(6u + Xn)b\ + 2X 11 (2n Q b\ + (alb\) + X ia a^ a - nb\ - 



2[iD h 
2^D\ n 



or equivalently, 



d 



XoiPoa^a - 
XoiOioa'a - 



where A is defined in Eq. I|34|) . To calculate (b + b^) = /?o + ao 
= (db{/dt) in Eq. (ggj, to obtain 



(bi) 
(b\) 



(66) 
(67) 

(68) 

(bi + b\) in the steady state, we setting (dbi/dt) — 

(69) 



2rjB m 
2^4 



2[iD m 



XoiPoa^a 
Xoiaoa^a 



Then after a simple calculation, we obtain the measured mean signal 

V^<&! + b{) = [(« - i(Su + A n + 4An|/3 | 2 ) + 2 J A 11 (a ) 2 ) #,] - h.c.} (o+o) 

= -iv 7 ^^- [K(/3o-a )-i(^ + A 11 + 2A 11 |/3 | 2 )(/3 + a )] (o+o). 



A 2 

Using Eq. (|32(l . we can simplify Eq. (|70() further and obtain 



y % { b 1 +b\)=-^2jl^(a^a). 



(70) 



(71) 



We note that the coefficient on the right hand side of Eq. (|71|l is proportional to VT, with a proportionality factor given 
by —y/d>n/ k{2N\ + 1). As the actual readout current is simply proportional to the average position of the ancilla 
oscillator 0, Eci. (|?Tj) gives the expected proportionality between the average measured current and the average phonon 
number of the system oscillator. 

In a typical experimental run, the measured current will contain a noise component made up of thermo-electrical 
noise in the transducer circuit as well as intrinsic quantum noise that arises directly from the back action noise 
when we measure phonon number. In order for the measurement to be quantum limited, we need to ensure that the 
dominant source of noise is back action noise. Recently, considerable progress towards this limit has been made in a 
nanoelectromechanical system l2lt 



VII. CONCLUSIONS 

We have investigated a scheme for the QND measurement of phonon number (cf ||) using two anharmonically 
coupled modes of oscillation of mesoscopic elastic bridges. We have included the self-anharmonic terms neglected in 
the previous analysis [{|, and analyzed the effect of higher order anharmonic terms in the approximation that the 
ancilla oscillator is heavily damped. We have shown that in the presence of self-anharmonic term, xf, of the ancilla 
oscillator, the effect of increasing driving strength and sclf-nonlinearity tends to shift the resonance frequency, increase 
the peak value and decrease the width of the response of the peak of (r/r ) as shown in Fig- El If the damping of 
the ancilla oscillator is much larger than the effect of the self-anharmonic term, the overall effect of self-anharmonic 
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term on the phonon number measurement is small for small detuning, justifying the assumption of neglecting the 
sclf-anharmonic term at zero detuning in Ref. [5j. Our calculation, however, allows one to do a quantitative analysis 
at finite detuning and without making this assumption. 

The key idea of the measurement scheme is that, from the point of view of the ancilla oscillator, the interaction 
with the system oscillator constitutes a shift in resonance frequency that is proportional to the time-averaged phonon 
number or energy excitation of the system oscillator. This frequency shift may be detected through a phase sensitive 
readout of the position of the driven readout oscillator. In a magnetic field, a wire patterned on the moving readout 
oscillator will result in an induced current which can be directly monitored by electrical means |8J. The current 
gives direct access to the position of the ancilla oscillator and, through the mechanism described in this paper, to 
the phonon number of the measured system oscillator, even in the presence of the self-anharmonic terms. We have 
shown that this scheme realizes an ideal QND measurement of phonon number in the limit that the back action 
induced phase diffusion rate is much larger than the rate at which transitions occur between phonon number states, 
(r/i^) — > oo. When the ratio (T/v) is finite and large, it is then possible to observe, in the readout current, quantum 
jumps between Fock (number states) in a mesoscopic mechanical oscillator, as the mechanical oscillator exchanges 
quanta with the environment. 

We briefly discuss below some possible realistic values for T and v. The value of Tq depends on external driving, 
as well as materials and dimensions of the mechanical beams (oscillators). Here we quote the example in Ref. [5j 
using two GaAs mechanical oscillators with resonance frequencies ljq = 2.3 GHz, u)\ = 0.36 GHz, and Q-factors 
Qo = 10000, Qi = 1000. The dimensions of the system oscillator are 0.6 /jm x 0.04 /im x 0.07 /j,m and those of the 
ancilla oscillator are 0.6 /im x 0.04 /j,m x 0.01 ^m. With the magnetic field 10 Tesla and the driving current 1 /lA, 
T and v will be T ~ 1.5 x 10 4 /s and v sa 1.2 x 10 6 /s, or T^/v — 0.013. A clear observation of quantum jumps 
requires Tq/i> 1, so that the present example is two orders of magnitude below the desired parameter regime. To 
increase the ratio of T to v we can improve on some of the parameters. One way is to increase the Q-factor of the 
system oscillator. Another way is to use lower density material such as carbon nanotubes as well as to decrease the 
thickness of the oscillator. These improvements are feasible with current fabrication technology. In addition, it is also 
possible to engineer the nonlinear coupling between the oscillators Q. Furthermore, different driving and detection 
schemes other than magnetomotive detection can be considered to increase the driving strength. Given the steady 
improvement in the fabrication technology and experimental techniques, we believe that observing quantum jumps 
between phonon number states in a mesoscopic oscillator will be possible in the near future. 
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APPENDIX A: EXPRESSIONS FOR THE TWO-TIME CORRELATION FUNCTIONS 

The two-time correlation functions in the main text for C (t, t'), where t > t' , are 

(ft (t) ft (<')) = 2^ {( a i + c ) cx p [~ A + (* - 0] + ( A i - c ) cx p M- (* - *')]} (ft 2 ) 

+ ^ {exp [-A+ (t - t')] - exp [-X- (t - t')]} (ftai) , (AI) 

( ai (t) ft (f )) - JL {(Ai + c) exp [-A+ (t - t')} + (A, - c) exp [-A_ (t - f )]} («iA> 

+ ^ {exp [-A+ (t ~ t')} - exp [-A_ (t - t')}} (af) , (A2) 

(ft (t) at (t')> = {exp [-A+ (t - t')} - cxp [-A_ (t - t')}} (ft 2 ) 

+ 7TT- {(Ai - c) exp [-A+ (i - t')} + (Ai + c) exp [-A_ (t - t')]} (ftai) , (A3) 
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\ 2 

(at (t) a x [i')) = {exp [-A+ (t - t')] - exp [-A_ (t - t')}} (aifr) 

+ ^ {(Ai - c) exp [-A+ (t - t')} + (Ai + c) exp [-A- (i - i')]} (a?) , (A4) 

where c = 4Anno +<5w + An and = (8lo + An) 2 + 8 (5ui + An) Anno + 12Af 1 nQ as in the main text. These equations 
give the c-number two-time correlation functions we need to obtain the operator two-time correlation functions in 
Eq. g2Hinj- 
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